Band structures of bilayer graphene superlattices 
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We formulate a low energy effective Hamiltonian to study superlattices in bilayer graphene (BLG) using a 
minimal model which supports quadratic band touching points. We show that a one dimensional (ID) periodic 
modulation of the chemical potential or the electric field perpendicular to the layers leads to the generation of 
zero-energy anisotropic massless Dirac fermions and finite energy Dirac points with tunable velocities. The 
electric field superlattice maps onto a coupled chain model comprised of 'topological' edge modes. 2D su- 
perlattice modulations are shown to lead to gaps on the mini-Brillouin zone boundary but do not, for certain 
symmetries, gap out the quadratic band touching point. Such potential variations, induced by impurities and 
rippling in biased BLG, could lead to subgap modes which are argued to be relevant to understanding transport 
measurements. 



Superlattices provide a route to band structure engineering 
in semiconductors Ut]. In graphene fl, a superlattice (SL) 
potential has been shown to lead to anisotropic Fermi veloc- 
ity renormalization [3], and generation of new Dirac points in 
the spectrum lHI-Ql resulting from the chiral nature of mass- 
less Dirac excitations. Such graphene SLs have been studied 
by epitaxial growth of graphene on Ir(lll) surface lUl Igl. Su- 
perlattice effects have also been studied in a topological in- 
sulator in proximity to a helical spin de nsity wave moH , and 
in graphene subject to a magnetic SL Ml ll 11211 . However, 
apart from transfer matrix studies of ID Kronig-Penney mod- 
els iHIlltl, SLs in bilayer graphene (BLG) have not been care- 
fully explored. 

Besides band structure engineering, there is a second moti- 
vation to study such BLG SLs. On theoretical grounds, BLG 
is an attractive candidate for transistor applications since it 
has a tunable gap which varies in proportion to the electric 
field perpendicular to the layers lll4[ llSll . However, trans- 
port measurements on BLG samples do not show the strong 
suppression of conductance at low temperatures expected on 
theoretical grounds lll4l lisil or from optical absorption mea- 
surements 1 16]. Instead, the transport data shows evidence for 
variable range hopping conduction [17-19] or a suppressed 
band gap il8ll20ll . It has been proposed that the observed ex- 
cess conductance arises from edge states [21], but transport 
measurements in a Corbino geometry do not support this sce- 
nario [j22(|, suggesting the existence of disorder-induced low 
energy modes in the bulk. To the extent that disorder poten- 
tials can be decomposed into Fourier components, we expect 
to learn something useful about disordered BLG by study- 
ing the simpler problem of periodic potential modulations in 
BLG. 

In this Letter, we study the band structures of BLG SLs, 
arising from periodic modulations of the chemical potential 
and the bias, using an effective low energy Hamiltonian. Our 
main results are the following, (i) Although the minimal 



model of BLG has quadratic band touching points, we find, 
remarkably, that a weak ID chemical potential modulation 
leads to the generation of linearly dispersing massless Dirac 
fermions with a tunable and anisotropic velocity. These Dirac 
fermion excitations are robust and rely on the chiral nature 
of the BLG quasiparticles. Beyond a critical modulation am- 
plitude, these Dirac modes get gapped out. (ii) An elec- 
tric field SL is shown to support linearly dispersing massless 
Dirac fermions and finite energy Dirac points which survive 
even for strong modulations. We provide a picture for these 
modes within a novel coupled chain model of 'topological' 
edge states, (iii) For 2D SLs, we show that for chemical poten- 
tial and electric field SLs the quadratic band touching points 
are protected for symmetric SLs with C4 or Cg symmetry, 
(iv) We compute the density of states for biased BLG with su- 
perimposed ID potential modulations, and find a plethora of 
subgap modes which we argue are important for understand- 
ing transport data. While our results on ID SLs overlap with 
work on Kronig-Penney models [13l 12311 . our analysis pro- 
vides simpler insights, highlights the role of the quasiparticle 
chirality, and is applied here to more general potential profiles 
as well as to 2D SLs. 

Effective Hamiltonian approach. — The low energy Hamil- 
tonian for Bernal-stacked BLG can be obtained by expand- 
ing its minimal tight binding spectrum near one of the Biil- 
louin zone corners (K points) lll411 . When the bias (i.e., inter- 
layer potential difference) is not too large, | A| <C t±,we find 
n = 4>^Hi^ [U], where 
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and ~ (flx, fex), with a (b) being the electron operator on 
the top (bottom) layer Here, tt = —id^ + dy, vp — \^td/2 ft! 
10^ m/s is the Fermi velocity, i?s3 eV is the nearest neighbor 
hopping integral, d sa 2.46 A is the distance between neigh- 
boring atoms on the same sublattice, Vi,2 are the potentials on 
each layer, and t± « 0.15t is the interlayer coupling. Unless 
stated, we set t=d^l. We will ignore inter-valley scattering 
assuming the potentials are varying slowly on the scale of d, 
so that identical physics is expected around the other valley (at 
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— K). Such an approach has been successfully used to study 
SLs in monolayer graphene f^^. 

To diagonalize i/kin^ we Fourier transform and then make 
a unitary transformation ap^{ap+|3p)/^/2, &p = e^'^p(ap — 

/3p)/ \/2, where cos 9p —Px/p and p = ^p'^ + p"^. This leads 

to 7?kin = Ep (ee(p)/3^/3p+e/i(p)a|[,ap)- Here ee,/i(p) = 
±p^/2m* are energies of electron (hole) states, with an ef- 
fective mass m* = t±/{2vp). This minimal model supports 
quadratic band touching points at ±K. 

When Vi^2(x) are periodic, we can also Fourier transform 
the SL potential to obtain iJsL = J2p g *^(p)W^p,g*(p - 
G), where 

^P'''-2\Vi{G)-V2{G)c^''> Vi{G) + V2{Gy J' 

vpt^p) ^ (q,!^^ )^ and 6 = 9p-G - dp is the angle between 
momenta p — G and p. Our aim is to understand the band 
structures of SLs described by ifkin + ^^sl- We will study 
ID SLs with period A along y, so that the reciprocal lattice 
vectors, {G}, are integer multiples of Q = (0,27r/A), and 
the mini Brillouin zone (MBZ) boundaries are at py = ±tt/X. 
We will also study 2D SLs. 

ID chemical potential superlattice. — Imposing a periodic 
potential 14(0;,?/) = V2{x,y) — C/(a;, y) corresponds to a 
chemical potential modulation. Numerically solving for the 
band structure of a periodic ID modulation using the above 
effective Hamiltonian, we find a pair of zero energy Dirac 
points in the MBZ in the vicinity of each valley. This is shown 
in FiglUfor a periodic step-like potential with (i) U {x, y) — U 
for < y < A/2 and (ii) U{x,y) = -U for X/2 < y < A. 
With increasing U, these Dirac points move away from each 
other along y. Beyond a critical modulation amplitude a full 
gap opens up. 

The existence of two Dirac cones at each valley is deeply 
rooted in the chiral nature of the low energy BLG quasiparti- 
cles, which causes the matrix elements of Eqn.|2]to depend on 
the scattering angle 6. For states with momenta parallel to the 
modulation direction, 9 — or ir, the off-diagonal matrix el- 
ements vanish; the electron and hole states then decouple, but 
electron-electron and hole-hole mixing is allowed. However, 




FIG. 1: Energy spectrum for a ID superlattice with step-like chemi- 
cal potential modulation of amplitude U. We set A = 60d, with [left 
panel] U — O.Olt showing two Dirac nodes split along y near K, 
and with [right panel] U = 0.04i: showing a full gap. 



in an extended zone scheme, all such electron (hole) states 
within the first MBZ only mix with electron (hole) states of 
higher (lower) energy, and so the energy of these states will 
be globally shifted down (up). This results in two level cross- 
ings along the modulation direction, which are protected by 
the chirality of the low energy BLG quasiparticles. If this 
electron-hole decoupling was true for all momenta, we would 
see the two parabolic bands crossing on a full circle in the 
MBZ, but going to momenta {6px,Py) leads to electron-hole 
mixing that is linear in Spx', this results in an avoided level 
crossing and the robust emergence of two Dirac cones in the 
MBZ. 

The location and velocity anisotropy of Dirac cones, as well 
as the critical modulation amplitude to gap them out, can 
be predicted using perturbation theory in U{G). The sec- 
ond order energy coiTection of states with p = (0,pj,) is 
A£;(2)(p) = En^o\UinCl)\V[eeMp)-eeAP + nQ)]. 
Since ee(p) < £e(p + nQ) while e/i(p) > e/i(p + nQ) in 
the MBZ, this correction is always negative (positive) for elec- 
tron (hole) states, as expected. 

Thus, the two bands will intersect and cross lin- 
early at momenta (0,±p*), where p*^/2m* = 

2™* T,n^o |t^("Q)lV ["^<3^ + 2p;rJ(3] . For weak modu- 
lations, p'y/Q 1, and keeping only n — ±1, we estimate 

p*y w V2m*|C/(Q)|A/7r. For a step profile, |C/(Q)| = 2U/tt, 
and |n| > 1 contributions are small. 

For small Spx away from the level crossing point, we can 
estimate the electron-hole mixing term using perturbation the- 
ory (|2J], and we find that the resulting eigenstates have en- 
ergies Cp = ±{16m*\UiQ)\^/\Q\^)Spx/p*y. The crossing 
points at (0, ±Py) are thus really massless Dirac points in the 
full MBZ. We find velocities Vy = p*y/m* w ^/2A|C/(Q)|/7r, 
and Vx = 2vy for the anisotropic linear dispersion. 

Once these Dirac nodes reach the MBZ boundary, Bragg 
scattering between them opens up a full gap. The critical po- 
tential strength, |[/c(Q)| for this is roughly estimated by set- 
ting p* = Q/2, which yields |C/c(Q)| ~ tt"^ /{y/2m*\^). For 
a step profile, with A = 60(i, we find Uc ~ 0.03i which is 
close to the numerical result 0.02t. 

ID electric field superlattice. — An electric field SL coiTe- 
sponds to Vi(a;, y) = —V2{x,y) = U{x,y). Solving for the 
resulting band structure, we find that it depends sensitively on 
the modulation type. To illustrate this, we consider a periodic 
potential, with U{y) = 2U{1 — w/X) for < y < w, and 
U{y) = —2Uw/X for w < y < X. We have set the average 
potential on each layer to be zero. If w = A/2, the resulting 
symmetric SL is found to support a pair of anisotropically dis- 
persing massless Dirac fermions at zero energy at (±p* , 0), as 
seen in Fig.|2](left panel). In addition, as shown in Fig.|2](iight 
panel), it supports a Dirac point at nonzero positive (as well 
as negative) energies at (0, tt/A) (or equivalently (0, — tt/A)). 
However, an asymmetric SL, with w ^ X/2, leads to a gap 
for all these Dirac fermions. More generally, we find that if 
the SL potential commutes with a generalized parity operator, 
V, which corresponds to y — y followed by exchanging the 
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FIG. 2: Energy spectram for a ID symmetric (see text) electric field 
superlattice with A = 60d and U = O.OSt, showing a pair of zero en- 
ergy massless Dirac fermions at (±p* , 0) [left panel] and a nonzero 
energy Dirac point at (0, ±7r/A) [right panel]. 



two layers of BLG, then these gapless Dirac points survive. 
Breaking V leads to gaps. 

A simple route to understanding these results that leads to 
other interesting predictions is to view the SL as a periodic 
array of 'kinks' and 'antikinks' where a kink (antikink) cor- 
responds to where the electric field flips from pointing up 
(down) to pointing down (up). A single such kink/antikink 
in the bias is well understood 121, .25n27il . In the absence of 
interactions a kink (antikink) supports a pair of light-moving 
(left-moving) 'topological' edge states near the K point for 
each spin. By time-reversal, these right and left movers get 
interchanged at the — K point. These modes are depicted in 
Fig. [3] (Although these modes were suggested to be topologi- 
cally protected, they are not truly stable against disorder; nev- 
ertheless disorder induced backscattering is weak li2l]| .) At a 
kink, we denote the higher (lower) energy edge state as tt (0), 
while we denote these states as tt (0) at an antikink. Hence, 
there are four points at each valley where kink and antikink 
modes cross: two of these occur at zero energy (tt-O and tt-O 
crossings), and two of them occur at nonzero energy (tt-tt and 
0-0 crossings). We will show below that these crossing points 
evolve into massless Dirac fermion modes in the MBZ of the 
SL. In order to see this, we construct a tight-binding model of 
such coupled 'topological' edge states. 

We observe that the Hamiltonian with the single kink (or 
antikink) potential is invaiiant under V, since H{y)V = 
<JxH{~y)(Jx — H{y). The 0/0 states are even under T', while 
the tt/tt states are odd under V ll25ll . Let us then construct 
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FIG. 3: (color online) Left: Spectrum of isolated kink (thin, red) and 
antikink (thick, blue). Higher (lower) energy modes are labelled vr 
(0) at a kink and as tt (0) at an antikink. Right: Schematic of hopping 
between the vr — vf and () — vr states. 



a reduced Hamiltonian which describes the hybridization be- 
tween neighboring edge modes. 

We begin with neighboring tt-O modes at zero energy and 
at a momentum p* (away from K). The hopping between 
neighboring 'wires' along y is then between states which have 
opposite velocities (since it is between a kink and an antikink 
edge state) and it is between a p-wave like state {V-odd) and 
an s-wave like state (T'-even). Using the index n to label the 
wires, the interchain hopping parameter will then alternate as 
{—l)^g for equally spaced wires and as q+5, —q+5 (with 5 < 
g) if pairs of wires are closer to each other |!24l|. Linearizing 
the dispersion at the crossing point, and letting vq denote the 
velocity of the linearized modes, 

n 

- J2(gi~ir + S){cl^Cp^n+i + h.c.) (3) 

n 

where p* is the location of the tt — crossing point 
in the single kink or antikink problem, and Cp^„ an- 
nihilates an electron on wire n with momentum p^- 
Let ^{px) = i^o{p ~ P*x)- Fourier transforming. 



we 



find 



HiPx) = Ep, ■ where 



HPx) = iC{px),-2gsin{py),-26cos{py)), with 



-Py 



^tt)'^, and Y^'p runs over the IVIBZ. The dispersion is 



thus E = ±^^'^{px) +452 cos'^(py) +4^2 sin^ipy). Conse- 
quently, when w = A/2, and the Hamiltonian commutes with 
V, we have (5 = and a Dirac cone is generated at (p* , 0), 
consistent with numerical results. When w 7^ A/2, the Hamil- 
tonian breaks V — we then have S 0, which leads to a gap 
4(5. Similar arguments hold for the other zero energy band 
crossing points. The velocity of the Dirac fermions is highly 
anisotropic and depends on g — this can be controlled by tun- 
ing the SL period and amplitude. 

The above analysis can also be repeated for the nonzero 
energy (0-0 and tt-tt) crossings |24J; in the symmetric case, 
w — A/2, we find Dirac cones at (0,±tt/A) on the MBZ. 
Once again, a modulation with w ^ X/2 results in a finite 6 
and opening of band gap. 

Interestingly, just as in poly acetylene, a domain wall be- 
tween a gapped region with w > X/2 and a gapped region 
with w < X/2 leads to new subgap soliton modes. Since 
each kink/antikink is itself like a domain wall, these should be 
viewed as soli tons in a soliton lattice! 

2D superlattices. — We have also considered 2D chessboard 
like SLs with fourfold rotation symmetry. For both types of 
2D SLs, chemical potential or electric field, the quadratic band 
touching point remains intact when the SL potential is 'sym- 
metric', ^1^2(2^ + A/2, y) = 1/1,2(2:, J/ + A/2) = ~Vi^2{x,y). 
This is consistent with the fact that no Dirac points can be 
generated in a way that conserves both topological charge and 
C4 (or Ce) symmetry 1281]. For asymmetric SLs, higher order 
corrections lead to modifications to the energy spectrum at the 
K-point [24]. For chemical potential SL, the charge neutrality 
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FIG. 4: (color online) Density of states for BLG subject to a uniform 
bias of A — O.lt and various chemical potential (left) and electric 
field (right) superlattices with period A = 60d. 



point (CNP) shifts slightly in energy, due to higher order ef- 
fects which reflect particle-hole symmetry breaking. For elec- 
tric field SLs, breaking generalized parity opens a small gap 
at the K-point [24.1 . 

Experimental implications. — Our work demonstrates that 
SL modulations in BLG can generate new Dirac fermion 
modes. Such modes are perturbatively stable to interaction ef- 
fects, and could be experimentally explored by suitable choice 
of substrates. Disorder will also lead to such bias and chem- 
ical potential modulations, albeit in random fashion. One 
source of such fluctuations is the presence of charged impu- 
rities, embedded in the underlying substrate (Si02) or, in the 
case of suspended BLG, in the residue of the etching/washing 
process. Such impurities are expected to locally shift the CNP, 
and to suppress or enhance the bandgap depending on the rel- 
ative sign of the bias and the impurity electric field [29]. If 
the impurity lies close to the surface it can locally reverse the 
parity of the interlayer bias leading to 'topological' subgap 
modes. Another source of SL fluctuations is rippling [30, 3ll, 
which would modulate the electric field perpendicular to the 
bilayer at the ripple wavelength. 

As a starting point to understanding the expected role of 
chemical potential and electric field fluctuations, Fig.|4]shows 
density of states (DOS) plots of a biased SL with periodic 
ID modulations. In the absence of a SL, the DOS diverges 
as 1/ ^/E at the gap edge arising from the ^ p"^ dispersion of 
modes near the gap edge. We find that both chemical potential 
or bias modulations, cause low energy subgap modes states in 
this system that will renormalize the average band gap, con- 
sistent with experiment. For chemical potential modulations, 
the subgap states are due to the local shift in the CNP. At fi- 
nite temperature, regions with a slightly shifted CNP will have 
thermally activated 'electron-hole' puddles that contribute to 
transport. For bias modulations, weak modulations locally 
enhance or suppress the bandgap, while strong modulations 
form 'topological' states in the bulk along interfaces where the 
field reverses sign [21 , 25-27]. The energy of these 'topologi- 
cal' midgap states decreases for large and dilute fluctuations, 
as the overlap between edge mode wavefunctions is reduced. 

Random potential fluctuations will have two important ef- 
fects not captured in our study of periodic modulations. First, 
it will cause the low energy density of states to broaden, caus- 



ing further suppression of the bandgap predicted by the peri- 
odic modulation. Second, dilute localized 'topological' states 
induced in the bulk by strong random electric field modula- 
tions due to charged impurities will contribute to transport 
through variable range hopping — this is broadly consistent 
with the tem pera ture dependence of the resistance in biased 
BLG lEMllll. 
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Note Added: After submission of this Letter, we received 
a preprint of Ref. ll32ll . which studies Dirac fermions in ID 
chemical potential superlattices in BLG and contains results 
consistent with ours. 
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ID CHEMICAL POTENTIAL SUPERLATTICE: 
PERTURBATION THEORY 



perturbative result is not strictly valid in this regime, it 
provides a crude estimate of 



Location of the Dirac Point — For small bias strengths, 
the location of Dirac point, (0,p*), can be estimated by 
setting the absolute value of the second order energy cor- 
rection equal to the free electron energy, 



1 



+ 



2Q2 _ 2p*ynQ 



(1) 



Equation [T] can then be expanded in the Py/Q ^ 1 limit 
to give 



(2) 



or retaining only the fundamental harmonic 



Pv 



V2m*\U{Q)\X 



(3) 



For the case of a square potential where |?7(nQ)| ~ 
2U/mr for odd n and keeping only the fundamental har- 
monic. 



Pv 



2V2m*UX/ 



(4) 



Alternatively, Eq. [2] can be computed explicitly through 
the relation V"'' = 7^4/96 to give pi « ± v^Ii|l£A. 

A comparison of the two expressions confirms that higher 
order harmonic corrections are indeed small (~ 1%) and 
can be ignored. 

Critical Superlattice Strength — Using the approxima- 
tion provided for p* in Eq. |3l it is also possible to make a 
crude estimate of the critical potential, Uc, before the on- 
set of a bandgap. From numerical results consistent with 
the perturbative results above, the Dirac points are seen 
to move out towards the mini-zone boundary with in- 
creasing U. Upon reaching the mini-zone boundary, any 
further increase in U results in an opening of a bandgap 
at the mini-zone boundary. Hence, the critical potential 
strength, Uc, can be determined from the condition that 
the Dirac point p*y — Pc ~ ±Q/2. Although the above 



\Ucm 



(5) 



V2m*A2' 

In the case of a step potential with A = 60fi, Uc ~ 

3 

— -tS- — T ~ O.OSi, which is reasonably close to the ob- 
served numerical value of 0.02i. 

Velocity Anisotropy — Also of interest is the degree of 
anisotropy of the group velocity about the Dirac cone. 
Along the p = iO,Py) direction, the above perturbation 
theory indicates that Vy = p*y/m* sa v^lM^Mi. To cal- 
culate the velocity along the Px direction, we perform a 
degenerate perturbation theory for an electron and hole 
states with momentum p = {pl0,pt)) for a small angle 
9 while retaining only the leading order harmonic. For 
such states, the off-diagonal terms in the Hamiltonian 
reduce to, 



p,p Q 



UiQ) 



*Msin(20p,p_Q) 



MQ)sin(20p,p_Q)i 

UiQ) r 



Given the small finite momentum along the px direc- 
tion, the originally pure electron states get mixed with 



and. 



l*p> = l/3p) + 



vice versa. To leadin 


g order in 9, 




c^(Q)l"p+Q> 


C/*(Q)|«p- 


Q> 


ee(p)-ee(p + Q) 


feip) - ee(p- 


Q) 


z0{/(Q)|/3p+Q> 


t9U*{Q)\/3p 


-q) 


ee(p) - e,,(p + Q) 


ee(p) - e/i(P 


Q) 


C^(Q)l/5p+Q) 1 




q) 


e/i(p) -e/i(P + Q) 


e/i(p) - e/i(p 


Q) 




i9U*{Q)\ap 


-Q> 


e/i(p) - ee(p + Q) 


eh(p) - ee(p 


-Q) 



to first order in 



Since, (vJ-^IFlM/p 

9 and are degenerate, we must compute the off-diagonal 
matrix elements {'iip\H\'^p) of the Hamiltonian. To lead- 
ing order in 9, 



{n\H\<)^-2mm' 



1 



1 



e(p) -e(p + Q) 
1 



6(p)-e(p-Q) e(p)-He(p + Q) 
1 



e(p) + e(p-Q; 



(7) 



2 



where e(p) = ee(p) = —^h{p)- In the small p limit, the 
matrix element reduces to 



mm*0\U{Q)\'^ 



(8) 



This give the solution Cp = ±16to*6I|{7(Q)P/|Q|2 to the 

perturbed Haniiltonian is then, from which the velocity 
can be calculated. Using the relation p*0=p^ for small 
the theta, 



e(Px) = 



16m*\U{Q)\^p, 
p;|Q|2 

= 2\/2A|C/(Q)|/ TT 



(9) 



Hence, the anisotropy of the velocity in the Dirac cone 
is predicted to be v^/vy = 2 for small U, which is again 
remarkably consistent with the numerical results. 

The above results are based on the two-band reduced 
Hamiltonian for BLG. We have also carried out a similar 
calculation for the four-band Hamiltonian. We found the 
same behavior when superlattice potential is not very 
strong. Surprisingly, when the superlattice potential is 
comparable to t±, more band touching points emerge in 
the mini Brillouin zone. 



ID ELECTRIC FIELD SUPERLATTICE: 
EFFECTIVE MODEL 

In this section, we demonstrate in detail how the Dirac 
cones generated by a symmetric bandgap modulation 
{w = \/2) can be understood in terms of a tight-binding 
theory that describes a chain of coupled ID wires. 

Transfer Integrals — To derive the relation between 
the transfer integrals connecting wire n and n + 1, and 
n — 1 and n, careful consideration of the symmetry of 
the soliton wavefunctions must be made. As stated in 
the main text, the principle symmetry of the wavefunc- 
tions follow from the invariance of the Hamiltonian un- 
der the combined operation of layer inversion and re- 
flection about a kink (or antikink), i.e. 'P^H{y)V = 
axH{—y)ax = H{y). Solutions are then of the form 



( f{y) 
\9{y) 



fiy) 
f{-y) 



fiy) 
-fi-y) 



(10) 



with corresponding eigenvalues of +1 and —1 of the op- 
erator V, respectively. For the case of the zero energy 
band crossing points (between the tt- and 0-modes or the 
TT- and 0-modes at either K-point), the soliton wavefunc- 
tions of the kink wire have opposite V symmetry to the 
soliton wavefunctions of the two neighbouring anti-kink 
wires. In contrast, for the case of the finite band cross- 
ing points (between the tt- and 7f-modes or the 0- and 
0-modes at either K-point), the soliton wavefunctions of 
the kink wire have the sam,e V symmetry as the soliton 
wavefunctions localized to its two neighbouring anti-kink 



wires. As we will now show, the transfer integrals de- 
scribing the hopping between neighbouring wires along 
the array is dependent on whether the parity of the cou- 
pled modes is the same or opposite. 

For concreteness, let us consider the region where the 
TT- and 0-bands cross at the K point for a symmet- 
ric modulation with w = A/2. Let us set the y = 
point to be an anti-kink wire. The anti-kink 0-modes 
have wavefunctions of the form ^°(y) = {w{y),w{—y)y 
while its two neighbours' 7r-modes have wavefunctions 
of the form *'^(y + A/2) = {v{y + X/2),-v{-y - X/2)f 
and ^'^{y-X/2) = {v{y - A/2), -v{-y + \/2)f , respec- 
tively. For simplicity, we will assume that the wavcfunc- 
tion overlap is finite only in the region between the wires, 
as this assumption does not effect our main result. The 
transfer matrix that determines the hopping parameter 
between the central wire and its left neighbour is then 

91= f ^°\y)H{y)W{y + X/2)dy, (11) 

J-A/2 

and between its right neighbour 



92 



A/2 



^''\y-\/2)H{y)^~\y)dy. (12) 



Inserting H{y) = axH{-y)ax and changing y ^ -y 
in the expression for g2, gives us the relation 



91 = -9*2 = \9W' 



(13) 



between alternate bonds. The same relationship holds 
for the other zero energy band crossing point. However, 
repeating this calculation for the finite energy band cross- 
ing points (where the wavefunctions have the same par- 
ity) yields the corresponding relation 



g'l 



92 



t : 



\9V 



(14) 



Without loss of generality, in both cases the hopping 
parameter be assumed to be real, as it is always possible 
to remove the phase factor by a simple gauge transforma- 
tion. Hence, we have deduced from very general symme- 
try arguments that when w = X/2 the hopping parameter 
between the 0- (0-) and tt- (tt-) modes of neighbouring 
wires alternates sign along the chain, and is uniform be- 
tween 0- (tt-) and 0- (tt-) modes of neighbouring wires. 

If we generalize this case where w ^ X/2, the sep- 
aration between neighboiiring wires is unequal and the 
magnitude of the hopping parameter will begin to alter- 
nate along the bonds. Again considering the tt- ()-band 
crossing at the K point and taking the the wire at larger 
y to be further than to the anti-kink wire than the other 
neighbour, gi = —{—g + S) and g^ = —{g + S), where 
g > is the average magnitude of the hopping between 
neighbouring wires and (5 > is the deviation. Alter- 
natively, for the finite energy band crossing points, the 
hopping parameter can be shown to be = — {g — 6) 
and g2 = -{g + 5). 
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FIG. 1. (Color online) Variation of the hopping parameter 
along the wire array. Left: Hopping between zero energy 
modes of opposite parity. Right: Hopping between finite en- 
ergy modes with the same parity. Shape, orientation and sign 
of the wavefunctions are completely schematic and serve only 
to be illustrative of efi'ect of parity on the hopping integrals. 
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FIG. 2. (Color online) Energy spectra for symmetric 2D 
chemical potential (left) and electric field (right) superlat- 
tices. Here, A = 60d, w = 30d, U = 0.025t. 



Velocity Anisotropy — Along the wire direction, p = 
{px, 0), the group velocity is the same as that of a single 
kink anti-kink pair at the band crossing point. This ve- 
locity is only sensitive the details of the bias profile vifithin 
the unit cell, and independent of the modulation period. 
In contrast, along the modulation direction, p = (0,pj,), 
the velocity is given by the effective interwire hopping 
strength, g, and is thus dependent on the period of the 
modulation, as g is given by the overlap between the 
wires. As a consequence, it is interesting that the veloc- 
ity along the each direction can be tuned independently. 
After tuning the velocity along the p = {px, 0) direction, 
the velocity along the p = {0,Py) direction can be tuned 
continuously by adjusting the period length while keep- 
ing the velocity along the p = {px, 0) direction fixed. 

Limits of Applicability — After having analyzed the 
low-energy model, we now comment on the limits in 
which this model can be applied. As with other tight- 
binding models, its validity is dependent on the extent 
of overlap between adjacent 'atomic-like' wavefunctions. 
In this system there are actually two ways the overlap 
can increase, either by decreasing the distance between 
adjacent wires or by extending the wavefunctions them- 
selves by decreasing the modulation amplitude. In the 
long period and/or large amplitude limit, the localized 
states decouple. This is akin to the atomic limit and 
explains the flattening of the dispersion along the Py di- 
rection. In addition, it is important to emphasize that 
in calculating the transfer integrals, we have made the 
approximation that the soliton wavefunctions are k- in- 
dependent. Hence, there is in fact some momentum de- 
pendence to the transfer integrals that has been ignored 
that may become relevant at higher energies and/or short 
modulation periods. 



2D SUPERLATTICES 

Band Structures of 2D BLG Superlattices — In this 
section, we will consider the band structure of BLG sub- 
ject to 2D square superlattices with period A in both 
directions. Specifically, we consider the chessboard like 



superlattice potentials, 

Vi{x,y) = Uu V2ix,y)^U2, 

< x,y < w, and w < x,y < X, 

Vi^2{x + X,y) = Vi^2{x,y + X) - 1^,2(2:, y), (15) 

and set the average of the superlattice potential to be 
zero in each supercell. For chemical potential super- 
lattices, Ui — U2 — U, and electric field superlattices, 
Ui = -U2 - U. 

Fig. [2] shows the highest valence band and the lowest 
conduction band of BLG subject to the above superlat- 
tices. Generally, the energy spectra are fully gapped on 
the MBZ boundaries. Around K point, the energy spec- 
tra depend on the details of the superlattices. When the 
following symmetry is present, 

Vi,2{x + A/2, y) = V,,2{x, y + A/2) = -^1,2(2:, y), (16) 

there is no gap opening at the K point. This can be 
understood from perturbation theory presented in next 
subsection. However, when w deviates from A/2, a gap 
will open at the K point for electric field superlattices. 
For chemical potential superlattices, no gap opens, but 
the charge neutrality point will shift in energy. 

Analysis of Spectrum at K Point — When superlat- 
tice potential Ui^2 is not very large, or the superlattice 
period A is not very small, we can understand the energy 
spectrum from perturbation theory. Let us focus on the 
electron states at if in a chemical potential superlattice. 
Up to second order, the energy correction from an elec- 
tron state with momentum G = (n, m) x 27r/A, where n 
and TO are integers, can be directly read off from the W- 
matrix, 



Ai?(2)(G) = 



1(1 



o2ieG 



)U{G) 



4(0-ee(G)) 



(17) 



where 9q, is the angle defined earlier. Similarly, the 
energy correction from a hole state with momentum 
G' = (— TO,ri) X 27r/A is 



(2). |(l-e2'»«')[/(G') 



(18) 
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Since Og, = 9g+tt/2, \U{G)\ = \U{G')\, and £e(G) = 
— £ft,(G'), the contribution from the above two states will 
cancel each other. Most importantly, due to the fourfold 
rotation symmetry, it is always possible to find such a 
pair of electron- hole states whose angles differ by 7r/2, 
which makes the cancellation exact even an infinite num- 
ber of states are taken into account. Therefore, the four- 
fold rotation symmetry dictates the second order energy 
correction to states at K to be zero, and this second 
order result does not depend on the details of the super- 
lattice. The above argument also applies for electric field 
superlattices. 

For symmetric superlattices, chemical potential or elec- 



tric field, there is no third order energy corrections from 
the superlattice potential due to the fact that U (G) = 
if either n or to is an even integer. Breaking of the sym- 
metry ()16|) . or more generally the symmetry associated 
with particle-hole transformation followed by translation 
by A/2, will lead to the modification of the energy spec- 
trum. Numerically, we have observed, in chemical poten- 
tial superlattices, the charge neutrality point will slightly 
shift in energy, a result we appear to recover in third or- 
der perturbation theory. In electric field superlattices, 
breaking of V is numerically found to lead to a gap at 
the K point. Further details of the band structure of 2D 
superlattices will be presented in a future publication. 



